###### Title ##########################
## Insecurity and Support for Female Leadership in Conflict States: 
## Evidence from Afghanistan by Jasmine Bhatia and
## Steve L. Monroe.
## Replication Command File: Conjoint Experiment (Corruption Placebo Analysis)  
## August 2023
##

# This file contains analysis and tables of the conjoint survey experiment
# presented in both the main article and supplementary information


### Table of Contents  ##################################
## 1. Load packages, data and clean data for MM analysis
## 2. Create and clean variables
## 3. H2: Analysis with "dirty" control group
## 4. Placebo test with Corruption Primes
## 5. Mechanisms: Insecurity and support provincial gov (broad control group) 


## Note: All code run using R 4.0.3 GUI 1.73 Catalina build (7892)

######## 1. Load Necessary Packages and Data ##############################

library(cjoint)
library(dplyr)
library(tidyr)
library(xtable)
library(tidyverse)
library(stats)
library(survey)
library(Rcpp)
library(estimatr)
library(knitr)
library(DT)
library(cli)
library(devtools)
devtools::install_github("m-freitag/cjpowR")
library(cjpowR)
library(stargazer)
library(MASS)
library(mlogit)
library("nnet")
library("scales")
library(here)

# set working directory and load dataset

getwd()


data <- read.csv("Afghan_Data_All_Clean.csv")


### 2. Format Data and Create and Clean Variables for Analysis #####

# Create separate lists of variables for each unique 
# job profile 1-8

profile_1 <- c("Rd_1_A_Gender", "Rd_1_A_Education", "Rd_1_A_Age", 
               "Rd_1_A_Ethnicity", "Rd_1_A_Professional_Experience",
               "Rd_1_A_Place_of_Birth")

profile_2 <- c("Rd_1_B_Gender", "Rd_1_B_Education", "Rd_1_B_Age", 
               "Rd_1_B_Ethnicity", "Rd_1_B_Professional_Experience",
               "Rd_1_B_Place_of_Birth")

profile_3 <- c("Rd_2_A_Gender", "Rd_2_A_Education", "Rd_2_A_Age", 
               "Rd_2_A_Ethnicity", "Rd_2_A_Professional_Experience",
               "Rd_2_A_Place_of_Birth")

profile_4 <- c("Rd_2_B_Gender", "Rd_2_B_Education", "Rd_2_B_Age", 
               "Rd_2_B_Ethnicity", "Rd_2_B_Professional_Experience",
               "Rd_2_B_Place_of_Birth")

profile_5 <- c("Rd_3_A_Gender", "Rd_3_A_Education", "Rd_3_A_Age", 
               "Rd_3_A_Ethnicity", "Rd_3_A_Professional_Experience",
               "Rd_3_A_Place_of_Birth")

profile_6 <- c("Rd_3_B_Gender", "Rd_3_B_Education", "Rd_3_B_Age", 
               "Rd_3_B_Ethnicity", "Rd_3_B_Professional_Experience",
               "Rd_3_B_Place_of_Birth")


# Create separate list of responses for each pair of profiles

responses_1 <- c("RD1_Con_Table_Choice", "RD_1_A_Rank", "RD_1_B_Rank")
responses_2 <- c("RD2_Con_Table_Choice", "RD_2_A_Rank", "RD_2_B_Rank")  
responses_3 <- c("RD3_Con_Table_Choice", "RD_3_A_Rank", "RD_3_B_Rank")    


# Subset relevant treatments, responses, and all covariates 
# for each profile

covars <- c("ControlGrp", "TreatInsc", "TreatBribe", "TreatNep",
            "Res_Gender", "Literacy", "HHhead",
            "Employ_Unem", "District",
            "Province", "Income", "Voted_Pres", "Res_Ethnicity",
            "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace") 

data2_1 <- data[, c("ID", "Rd_1_A_Gender", "Rd_1_A_Education", "Rd_1_A_Age", 
                    "Rd_1_A_Ethnicity", "Rd_1_A_Professional_Experience",
                    "Rd_1_A_Place_of_Birth",
                    "RD1_Con_Table_Choice", "RD_1_A_Rank", "RD_1_B_Rank",
                    "ControlGrp", "TreatInsc", "TreatBribe", "TreatNep",  "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace")]


data2_1$profile <- rep(1, 2485)                    


data2_1$Pair <- "A" #name of pair. Will use this to create a profile pair id

data2_2 <- data[, c("ID", "Rd_1_B_Gender", "Rd_1_B_Education", "Rd_1_B_Age", 
                    "Rd_1_B_Ethnicity", "Rd_1_B_Professional_Experience",
                    "Rd_1_B_Place_of_Birth",
                    "RD1_Con_Table_Choice", "RD_1_A_Rank", "RD_1_B_Rank",
                    "ControlGrp", "TreatInsc", "TreatBribe", "TreatNep", 
                    "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace")]


data2_2$profile <- rep(2, 2485)   

data2_2$Pair <- "A" #name of pair. Will use this to create a profile pair id

data2_3 <- data[, c("ID", "Rd_2_A_Gender", "Rd_2_A_Education", "Rd_2_A_Age", 
                    "Rd_2_A_Ethnicity", "Rd_2_A_Professional_Experience",
                    "Rd_2_A_Place_of_Birth",
                    "RD2_Con_Table_Choice", "RD_2_A_Rank", "RD_2_B_Rank",
                    "ControlGrp", "TreatInsc", "TreatBribe", "TreatNep", 
                    "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace")]


data2_3$profile <- rep(3, 2485)                    

data2_3$Pair <- "B" #name of pair. Will use this to create a profile pair id

data2_4 <- data[, c("ID", "Rd_2_B_Gender", "Rd_2_B_Education", "Rd_2_B_Age", 
                    "Rd_2_B_Ethnicity", "Rd_2_B_Professional_Experience",
                    "Rd_2_B_Place_of_Birth",
                    "RD2_Con_Table_Choice", "RD_2_A_Rank", "RD_2_B_Rank",
                    "ControlGrp", "TreatInsc", "TreatBribe", "TreatNep", 
                    "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace")]

data2_4$profile <- rep(4, 2485)  

data2_4$Pair <- "B" #name of pair. Will use this to create a profile pair id

data2_5 <- data[, c("ID", "Rd_3_A_Gender", "Rd_3_A_Education", "Rd_3_A_Age", 
                    "Rd_3_A_Ethnicity", "Rd_3_A_Professional_Experience",
                    "Rd_3_A_Place_of_Birth",
                    "RD3_Con_Table_Choice", "RD_3_A_Rank", "RD_3_B_Rank",
                    "ControlGrp", "TreatInsc", "TreatBribe", "TreatNep", 
                    "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace")]


data2_5$profile <- rep(5, 2485)                    

data2_5$Pair <- "C" #name of pair. Will use this to create a profile pair id

data2_6 <- data[, c("ID", "Rd_3_B_Gender", "Rd_3_B_Education", "Rd_3_B_Age", 
                    "Rd_3_B_Ethnicity", "Rd_3_B_Professional_Experience",
                    "Rd_3_B_Place_of_Birth",
                    "RD3_Con_Table_Choice", "RD_3_A_Rank", "RD_3_B_Rank",
                    "ControlGrp", "TreatInsc", "TreatBribe", "TreatNep", 
                    "Res_Gender", "Literacy", "HHhead",
                    "Employ_Unem", "District",
                    "Province", "Income", "Voted_Pres", "Res_Ethnicity",
                    "IntForces", "TrustProvGov", "TrustCentGov", "Behav_Peace")]


data2_6$profile <- rep(6, 2485)  

data2_6$Pair <- "C" #name of pair. Will use this to create a profile pair id

# Rename columns and rbind datasets back together
neutral_names <- c("ID", "Gender", "Education", "Age", 
                   "Ethnicity", "Career", "POB",
                   str_replace(responses_1, "_1", ""), 
                   covars, "profile", "Pair")

names(data2_1) <- neutral_names
names(data2_2) <- neutral_names
names(data2_3) <- neutral_names
names(data2_4) <- neutral_names
names(data2_5) <- neutral_names
names(data2_6) <- neutral_names


data2 <- data.frame(rbind(data2_1, data2_2, 
                          data2_3, data2_4, 
                          data2_5, data2_6 
))

# Add a unique row identifier
data2$num <- 1:14910

## Create Pair Variable

data2$PairID <- paste(data2$ID, data2$Pair, sep="")


## Put variables in correct format for cjoint package functions

# Rename data for easier typing: 
data <- data2 


# Consolidate rating variables
data$rating <- case_when(
                data$profile %in% c(1,3,5) ~ data$RD_A_Rank,
                data$profile %in% c(2,4,6) ~ data$RD_B_Rank
)


# Make a dummy variable for choice  
# Note: NAs must be coded as zero for cjoint functions to work

data$choice_dummy <- ifelse(
                      data$profile %in% c(1,3,5) & data$RD1_Con_Table_Choice == 1 |
                      data$profile %in% c(2,4,6) & data$RD1_Con_Table_Choice == 2, 1, 0)


# Military Variable


data$Military <- ifelse(data$Career == "Military", "Military", "Civilian")
data$Military <- as.factor(data$Military)


# Gender Attribute Variable

data$Gender <-as.factor(data$Gender)

# Outcome Varibles

data$Choice <- data$choice_dummy
data$Rating <- data$rating


# Create Treatment Variable
# Make Variable a Dummy: Insecurity vs. Everything Else

data$TreatGroup <- data$TreatInsc

# Set NAs of TreatGroup to 0

data$TreatGroup[is.na(data$TreatGroup)] <- 0

data$TreatGroup <- ifelse(data$TreatGroup > 0, 1, 0)

data$TreatGroup2 <- ifelse(data$TreatGroup == 1, "Insecurity Prime", "Other Text")

data$TreatGroup2 <- factor(data$TreatGroup2, 
                     levels = c("Other Text", "Insecurity Prime"))

## Make factorial variable for the different treatments

# Create variables for each

data$TreatInsc1 <- data$TreatInsc

data$TreatInsc1[is.na(data$TreatInsc1)] <- 0

data$TreatInsc1 <- ifelse(data$TreatInsc1 > 0, 1, 0)

data$TreatNep1 <- data$TreatNep

data$TreatNep1[is.na(data$TreatNep1)] <- 0

data$TreatNep1 <- ifelse(data$TreatNep1 > 0, 1, 0)

data$TreatBribe1 <- data$TreatBribe

data$TreatBribe1[is.na(data$TreatBribe1)] <- 0

data$TreatBribe1 <- ifelse(data$TreatBribe1 > 0, 1, 0)

data$Control1 <- data$ControlGrp

data$Control1[is.na(data$Control1)] <- 0

data$Control1 <- ifelse(data$Control1 > 0, 1, 0)


data$Treat_Factor <- ifelse(data$TreatInsc1 == 1, "Insecurity",
                     ifelse(data$TreatNep1 == 1, "Nepotism",      
                     ifelse(data$TreatBribe1 == 1, "Bribery",
                      ifelse(data$Control1 == 1, "Control",
                         NA))))      


# Make a Variable for Treatment and Gender

data$Respondent_Gender <- ifelse(data$Res_Gender == 1, "Male",
                                 "Female")

data$TreatGender <- ifelse(data$TreatGroup == 1 & data$Respondent_Gender == "Female",
                           "Insecurity Prime - Female",
                           ifelse(data$TreatGroup == 1 & data$Respondent_Gender == "Male",
                           "Insecurity Prime - Male",     
                             ifelse(data$TreatGroup == 0 & data$Respondent_Gender == "Male",
                              "All Other Text - Male", "All Other Text - Female")))       


data$TreatGender <- as.factor(data$TreatGender)


# Treat Gender 3 (Merge the corruption variables )

data$TreatGender3 <- ifelse(data$Treat_Factor == "Insecurity" & data$Respondent_Gender == "Female",
                       "Insecurity Prime - Female",
                       ifelse(data$Treat_Factor == "Insecurity" & data$Respondent_Gender == "Male",
                         "Insecurity Prime - Male",     
                           ifelse(data$Treat_Factor == "Nepotism" & data$Respondent_Gender == "Female",
                            "Corruption Prime - Female",
                             ifelse(data$Treat_Factor == "Nepotism" & data$Respondent_Gender == "Male",
                            "Corruption Prime - Male",
                            ifelse(data$Treat_Factor == "Bribery" & data$Respondent_Gender == "Female",
                             "Corruption Prime - Female",
                             ifelse(data$Treat_Factor == "Bribery" & data$Respondent_Gender == "Male",
                            "Corruption Prime - Male",
                           ifelse(data$Treat_Factor == "Control" & data$Respondent_Gender == "Male",
                          "Neutral Text - Male", "Neutral Text - Female")))))))

data$TreatGender3 <- as.factor(data$TreatGender3)



# Make a new international forces support variable

data$IntForces_new <- case_when(
  data$IntForces == 1 ~ 5,
  data$IntForces  == 2 ~ 4, 
  data$IntForces == 3 ~ 3, 
  data$IntForces == 4 ~ 2,
  data$IntForces == 5 ~ 1,
  TRUE ~ NA_real_
)


# Make support for central gov variable


data$TrustCentGov_new <- case_when(
  data$TrustCentGov == 1 ~ 5,
  data$TrustCentGov  == 2 ~ 4, 
  data$TrustCentGov == 3 ~ 3, 
  data$TrustCentGov == 4 ~ 2,
  data$TrustCentGov == 5 ~ 1,
  TRUE ~ NA_real_
)

# Make support for provincial government variable

data$TrustProvGov_new <- case_when(
  data$TrustProvGov == 1 ~ 5,
  data$TrustProvGov  == 2 ~ 4, 
  data$TrustProvGov == 3 ~ 3, 
  data$TrustProvGov == 4 ~ 2,
  data$TrustProvGov == 5 ~ 1,
  TRUE ~ NA_real_
)


# Make a new peace variable

data$Behav_Peace_new <- case_when(
  data$Behav_Peace == 1 ~ 6,
  data$Behav_Peace  == 2 ~ 5, 
  data$Behav_Peace == 3 ~ 4, 
  data$Behav_Peace == 4 ~ 3,
  data$Behav_Peace == 5 ~ 2,
  data$Behav_Peace== 6 ~ 1,
  TRUE ~ NA_real_
)



data$dove <- ifelse(data$IntForces_new < 3, "Dove", "Hawk")

data$dove <- as.factor(data$dove)

data$high_prov_trust <- ifelse(data$TrustProvGov_new > 3, 1, 0)

data$high_prov_trust_factor <- ifelse(data$high_prov_trust == 1, "High Trust",
                                      "Low Trust")
data$high_prov_trust_factor <- as.factor(data$high_prov_trust_factor)


data$high_center_trust <- ifelse(data$TrustCentGov_new > 2, 1, 0)

data$high_center_trust_factor <- ifelse(data$high_center_trust == 1, "High Trust",
                                      "Low Trust")
data$high_center_trust_factor <- as.factor(data$high_center_trust_factor)


data$high_peace <- ifelse(data$Behav_Peace_new == 6, "High Peace",
                          "Low Peace")

data$high_peace_factor <- as.factor(data$high_peace)



### 3. Analysis H2 with dirty control group (neutral text + corruption primes) ####


# Load cregg package fro MM analysis
# install.packages("cregg")
library("cregg") 
# note: cregg masks amce function from cjoint, 
# detach this package before attempting to 
# re-run amce analysis above or call cjoint::amce


## SI Table 10: Corruption Primes added to Control Group (SI Section 2.8)


# Calculate marginal means
h2_extended <- cj(data, Choice  ~ Gender,
                  id = ~ ID,
                  estimate = "mm",
                  by = ~ TreatGender
)


h2_extended <- h2_extended[, c("BY", "level", "estimate", "std.error", "TreatGender")]

estimate <- h2_extended[, "estimate"]

ses <- h2_extended[, "std.error"]

table <- rbind(estimate, ses)

## Female Preferences

table_female <- table[, c(1:2,5:6)]

table_female_1 <- table_female[, c(1,3)]
table_female_2 <- table_female[, c(2,4)]

table_female <- rbind(table_female_1, table_female_2)

h2_female <- round(table_female, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_female <- as.data.frame(cbind(Gender, h2_female))

colnames(h2_female) <- c("Leader's Gender", "All Other Text", "Insecurity Prime")

view(h2_female)


## Male Preferences

table_male <- table[, c(3:4,7:8)]

table_male_1 <- table_male[, c(1,3)]
table_male_2 <- table_male[, c(2,4)]

table_male <- rbind(table_male_1, table_male_2)

h2_male <- round(table_male, digits = 3) 

Gender <- c("Female", "", "Male", "")

h2_male <- as.data.frame(cbind(Gender, h2_male))

colnames(h2_male) <- c("Leader's Gender", "All Other Text", "Insecurity Prime")

view(h2_male)

# get observations

table(data$TreatGender)


# Hand code Table 10 into the SI

# Subset to Female Respondents

female <- data[data$Respondent_Gender == "Female",]


# Calculate marginal means
h2_extended_female <- cj(female, Choice  ~ Gender,
                         id = ~ ID,
                         estimate = "mm",
                         by = ~ TreatGroup2,
                         h = 0.5
)


## F-Test results reported in text
f.test.h2_extended_female <- cj_anova(female, Choice  ~ Gender,
                                      by = ~ TreatGroup2)

f.test.h2_extended_female # p value: 0.156 , reported in Manuscript Section 4.2



# Power Analysis: Male respondents with extended data


## Subset to Men only

male <- data[data$Respondent_Gender == "Male",]


amces_male <- cj(male, Choice ~ Gender, 
                 id = ~ID, estimate = "amce", 
                 by = ~TreatGroup2)


# amce for insecurity male: 0.081
# n = 6792

cjpowr_amcie(
  delta1 = 0,
  delta3 = 0.08,
  alpha = 0.05,
  n = 6792,
  levels1 = 2,
  levels2 = 2,
  p00= 0.375,
  p01=0.125,
  p10=0.375,
  p11=0.125,
  sims = 1e+05,
  delta0 = 0.5
)


# power: 0.82 (choice), Reported in manuscript section 4.2


## 4. Corruption Placebo and Plot (SI Section 2.12) #####


# Remove Insecurity Prime from Analysis 

data$Corrupt_Treat <- ifelse(data$TreatGender3 == "Insecurity Prime - Female", NA,
                       ifelse(data$TreatGender3 == "Insecurity Prime - Male", NA, 
                         data$TreatGender3))       


data$Treatment <-  ifelse(data$Corrupt_Treat == 5, "Neutral Text - Female",
                    ifelse(data$Corrupt_Treat == 1, "Corruption Prime - Female",
                    ifelse(data$Corrupt_Treat == 6, "Neutral Text - Male",     
                    ifelse(data$Corrupt_Treat == 2, "Corruption Prime - Male",      
                    NA))))

data$Treatment <- factor(data$Treatment,
                  levels = c("Neutral Text - Female",
                   "Corruption Prime - Female",
                    "Neutral Text - Male",
                     "Corruption Prime - Male"))


corruption_placebo <- cj(data, Choice ~ Gender, 
                         id = ~ID, estimate = "mm", by = ~ Treatment)


group <- c("Female Leader Neutral Text - Female Respondent", 
           "Male Leader Neutral Text - Female Respondent",
           "Female Leader Corruption Prime - Female Respondent",
           "Male Leader Corruption Prime - Female Respondent",
           "Female Leader Neutral Text - Male Respondent", 
           "Male Leader Neutral Text - Male Respondent",
           "Female Leader Corruption Prime - Male Respondent",
           "Male Leader Corruption Prime - Male Respondent")

corruption_placebo <- cbind(corruption_placebo, group)


corruption_placebo$group2 <- factor(corruption_placebo$group,
                    levels = c("Female Leader Corruption Prime - Female Respondent",
                               "Male Leader Corruption Prime - Female Respondent",
                               "Female Leader Neutral Text - Female Respondent", 
                               "Male Leader Neutral Text - Female Respondent",
                               "Female Leader Corruption Prime - Male Respondent",
                               "Male Leader Corruption Prime - Male Respondent",
                               "Female Leader Neutral Text - Male Respondent", 
                               "Male Leader Neutral Text - Male Respondent"
                                ))


# SI Figure 5 

si_figure5 <- ggplot(data = as.data.frame(corruption_placebo),
                   aes(x = estimate,
                   y = group2)) + 
                  geom_point(size = 3)+ 
                  geom_pointrange(aes(xmin= estimate - 2*std.error,
                  xmax=estimate + 2*std.error))+
                  labs(x="Marginal Mean", y = "")+
                  geom_vline(xintercept = 0.5, linetype = "dashed")+
                  scale_y_discrete(labels=c("Male Leader Neutral Text - Male Respondent" = 
                              "Male Leader\n(Neutral Text - \n Male Respondent)",
                            "Male Leader Corruption Prime - Male Respondent" = 
                              "Male Leader\n(Corruption Prime - \n Male Respondent)",
                            "Female Leader Neutral Text - Male Respondent" = 
                              "Female Leader\n(Neutral Text - \n Male Respondent)",
                            "Female Leader Corruption Prime - Male Respondent" = 
                              "Female Leader\n(Corruption Prime - \n Male Respondent)",
                            "Male Leader Neutral Text - Female Respondent" = 
                              "Male Leader\n(Neutral Text - \n Female Respondent)",
                            "Male Leader Corruption Prime - Female Respondent" = 
                              "Male Leader\n(Corruption Prime - \n Female Respondent)",
                            "Female Leader Neutral Text - Female Respondent" = 
                              "Female Leader\n(Neutral Text - \n Female Respondent)",
                            "Female Leader Corruption Prime - Female Respondent" = 
                              "Female Leader\n(Corruption Prime - \n Female Respondent)"))



# SI Table 14 Corruption Table

corruption_placebo_table <- corruption_placebo[, c("BY", "level", "estimate", "std.error", "Treatment")]

estimate <- corruption_placebo_table[, "estimate"]

ses <- corruption_placebo_table[, "std.error"]

table <- rbind(estimate, ses)

## Female Preferences

table_female <- table[, c(1:4)]

table_female_1 <- table_female[, c(1,3)]
table_female_2 <- table_female[, c(2,4)]

table_female <- rbind(table_female_1, table_female_2)

female_corruption <- round(table_female, digits = 3) 

Gender <- c("Female", "", "Male", "")

female_corruption <- as.data.frame(cbind(Gender, female_corruption))

colnames(female_corruption) <- c("Leader's Gender",  "Corruption Prime", "Neutral Text")

view(female_corruption)


## Male Preferences

table_male <- table[, c(5:8)]

table_male_1 <- table_male[, c(1,3)]
table_male_2 <- table_male[, c(2,4)]

table_male <- rbind(table_male_1, table_male_2)

male_corruption <- round(table_male, digits = 3) 

Gender <- c("Female", "", "Male", "")

male_corruption <- as.data.frame(cbind(Gender, male_corruption))

colnames(male_corruption) <- c("Leader's Gender", "Corruption Prime", "Neutral Text")

view(male_corruption)

# get observations

table(data$Treatment)


## SI Table 14, handcode into SI


### 5. Mechanisms: Insecurity and support provincial gov (broad control group) #####


# remove treated respondents
control <- data[data$TreatInsc1 == 0,]

control$dove <- as.factor(control$dove)


# subset to female respondents who did not receive insecurity prime

female_control <- control[control$Respondent_Gender == "Female",]


# Provincial support

control_women_prov <- cj(female_control, 
                           Choice ~ Gender, id = ~ID, estimate = "mm", 
                           by = ~ high_prov_trust_factor)


f.test.mechanism_prov_female <- cj_anova(female_control, Choice ~ Gender, 
                                         by = ~ high_prov_trust_factor)

## p: 0.002, f-test p value is < 0.01, reported in manuscript section 4.3


## SI Figure 1 in SI Section 2.2

group <- c("Female Leader (High Trust)", 
           "Male Leader (High Trust)",
           "Female Leader (Low Trust)",
           "Male Leader (Low Trust)")


h2_prov <- cbind(control_women_prov, group)

# Relevel the treatment  groups 
h2_prov$group2 <- factor(h2_prov$group,
                    levels = c("Female Leader (Low Trust)",
                               "Male Leader (Low Trust)",
                               "Female Leader (High Trust)",
                               "Male Leader (High Trust)"))


# Make SI Figure 1 

si_figure1 <- ggplot(data = as.data.frame(h2_prov),
                     aes(x = estimate,
                         y = group2)) + 
  geom_point(size = 3)+ 
  geom_pointrange(aes(xmin= estimate - 2*std.error,
                      xmax=estimate + 2*std.error))+
  labs(x="Marginal Mean", y = "")+
  #ggtitle("Trust in Provincial Governance and Support for \n Female Leadership (Women, Control)")+
  geom_vline(xintercept = 0.5, linetype = "dashed")



# Doves

h_control_women_dov <- cj(female_control, 
                          Choice ~ Gender, id = ~ID, estimate = "mm", 
                          by = ~ dove)


f.test.mechanism_dov_female <- cj_anova(female_control, Choice ~ Gender, 
                                        by = ~ dove)

## not stat for doves. 


## Insecurity and support for women among men with different
## levels of trust in provincial gov

male_control <- control[control$Respondent_Gender == "Male",]


# Provincial support

control_men_prov <- cj(male_control, 
                         Choice ~ Gender, id = ~ID, estimate = "mm", 
                         by = ~ high_prov_trust_factor)


f.test.mechanism_prov_male <- cj_anova(male_control, Choice ~ Gender, 
                                       by = ~ high_prov_trust_factor)

## F-test, p = 0.25, reported in Manuscript Section 4.3

